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PREFACE 


The  author  thanks  Robert  F.  Stellingwerf  for  introducing  him  to  the  SPH 
method  and  for  many  enlightening  discussions.  This  work  was  supported  under  Con¬ 
tract  DMA  00I-88-C-0079. 
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•The  becquercl  (Bq)  is  the  SI  unit  o(  radioactivity;  1  Bq  -  1  event/* 
••The  Gray  (G>>  is  the  SI  unit  of  absorbed  radiation. 
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SECTION  1 


INTRODUCTION 


Although  particle  methods  in  hydrodynamics  have  been  remarkably  suc¬ 
cessful  in  many  problems,'""4  current  methods  suffer  from  an  inadequate  treatment 
of  boundary  conditions.  This  is  particularly  evident  when  one  tries  to  include  heat 
transport  in  the  energy  equation.5"6  Without  accurate  values  of  the  flux  at  the  bound¬ 
aries,  it  is  not  possible  to  calculate  the  net  exchange  of  energy  between  a  fluid  and  its 
environ  rnent. 

Another  class  of  problems  requiring  accurate  treatment  of  boundary  condi¬ 
tions  is  illustrated  by  the  shock  tube  problem  where  the  shock  is  driven  by  a  piston. 

In  this  case  the  particle  method  must  be  able  to  accommodate  an  externally- 
applied  (boundary)  pressure.  The  method  should  also  be  able  to  handle  conduction 
heating  of  the  gas  by  a  hot  wall  or  piston;  and,  finally,  boundaries  must  be  treated 
accurately  to  include  radiative  heat  exchange,  in  particular,  radiative  cooling. 

In  this  paper,  some  new  algorithms  are  developed  for  smooth  particle  hy¬ 
drodynamics  in  problems  where  external  boundary  conditions  are  imposed.  Although 
thermal  and  radiative  diffusion  are  included,  radiation  transport  is  neglected  in  the 
present  paper. 
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SECTION  2 


KERNEL  ESTIMATES  NEAR  BOUNDARIES 


The  standard  kernel  estimate  of  a  function  is  given  by3 

f*{x)  -■  f  f{x')W(x  -  x',h)dx'  , 

J  0 

where  the  kernel  is  normalized  by 


In  these  expressions  the  domain  of  definition  is  0  <  x  <  X . 


As  long  as  h  — >  0  one  would  not  expect  normalization  problems  near  the 
boundary,  but  in  practical  applications  kernels  with  finite  h  are  used.  One  would  then 
expect  difficulties  within  x  ~  h  of  the  boundary  because  the  integrals  are  truncated  by 
the  boundary.  For  example 


lirn  J  W (x  -  x')dx'  =  ^  ,  (3) 


which  is  obvious  from  Fig.  1. 

Consider  what  this  does  to  the  kernel  estimate  of  the  density.  Assume  a  dense 
constant  spacing  of  particles  in  the  vicinity  of  the  boundary  X.  Far  from  the  wall,  the 
kernel  estimate  gives 


p{  x)  —  f  p{x')W  (x  —  x\h)dx 

■J  O 

~  m  W(x  —  Xj,h) 

; 

""  Pc  1  (-1) 


2 


I 


Figure  1.  A  typical  smoothing  kernel,  showing  how  the 
symmetrical  form  is  truncated  as  the  particle  ap¬ 
proaches  a  boundary. 

where  pc  is  the  constant  value.  However,  following  Eq.  3,  as  x  — ►  X ,  the  kernel  estimate; 
becomes 


P{X)  = 


lim  [  p{x)W [x  -  x  ,h)dx 

z^X  J0 

m  ^  W [X  -  Xj .  h) 


:Pc  ■ 


So  within  ~  /.  of  the  boundary  the  interpolation  feature  of  the  kernel  estimate  give's 
spurious  resuks.  Also,  presumably,  one  finds 


r(x) 


(fi) 


fo r  the  case  of  a  uniform  pressure  in  the  neighborhood  of  the  boundary  as  in  Fig.  2. 


3 


-Pc 
-1  12  P 

c 

h  X 


Figure  2.  The  falloff  in  a  uniform  pressure  distribution  as 
given  by  a  kernel  estimate  in  the  neighborhood 
of  the  boundary. 

In  summary,  we  have  seen  that  because  of  the  truncation  of  the  normalization 
integral  within  ~  h  of  a  boundary,  interpolation  implied  by  the  kernel  estimate  may 
not  be  accurate.  Renormalization  is  a  possibility  but  does  not  appear  to  work.  Also, 
the  boundary  value  method  presented  below  seems  to  require  the  apparent  failure 
manifested  in  Eq.  3. 
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SECTION  3 


MOMENTUM  EQUATION  WITH  BOUNDARY 

PRESSURE 


All  treatments  of  SPH  of  which  the  author  is  aware  require  integration  by 
parts  when  deriving  kernel  estimates  but  simply  drop  the  boundary  terms.  (See  Section 
8  below.)  This  is  presumably  because  the  kernel  is  supposed  to  mimic  a  (5-function 
which  goes  to  zero  sufficiently  far  from  the  particles  representing  the  edge  of  the  system. 

However,  if  we  are  to  consider  problems  where  the  particles  interact  with  a 
boundary  condition,  such  as  an  externally-applied  pressure,  we  must  specifically  include 
the  boundary  terms.  In  this  section  the  two  procedures  developed  by  Monaghan  for 
the  momentum  equation2,3  are  modified  to  include  boundary  pressure. 

In  the  first  case,  write  the  momentum  equation  as 

dv  1  dP 

dt  p  dx 

-4(3-3 

The  kernel  estimate  of  this  equation  is  formed  by  changing  the  independent  variable 
to  x',  multiplying  by  W[x  —  x')dx'  and  integrating  over  the  domain  0  <  x  <  A\  One 
obtains 


where  the  subscrip'  denotes  “kerne!  estimate.”  The  second  integral  has  been  lin¬ 
earized  by  evaUia’ing  (P  /  p‘)  by  its  kernel  estimate  and  taking  it  outside  the  integral 
sign.  This  is  justified  since  VV(x  --  x()  acts  like  6(x  —  x). 
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After  integrating  by  parts,  one  finds 


dvK 

dt 


-W{x 


'  '  X  rX  /  p\ 

l  o  \p) 


*  (  P\  dW{x  -  x') 


dx 


4)  [pW{x-x) | 
J  K 


i\  ]X 

O 


y)J> 


dx 

X  dw{ x-x') 


dx 


dx 1  , 


(9) 


where  we  have  used  dW/dx  =  —  dW/dx  assuming  a  symmetric  kernel. 

If  we  now  evaluate  the  integrals  by  the  particle  method,  replacing 


/ 


f(x')W(x  -  x')dx'  —>  -  x})  —  , 

i  P ; 


(10) 


we  obtain  an  expression  for  the  momentum  equation  with  external  boundary  pressures. 


dv{ 

dt 


P  j  (  P 
P\  IP 


L \p2J0  +  Vp2 


P2 )  x  +  \P2 


,P2)}+\P2 


p0W{xi  -  o) 

pxW(xt  -  X) 

dW{xt  -  Xj) 


dx, 


Thus,  particle  i  feels  the  effect  of  the  boundary  pressure  if  it  is  within  range  of  the 
boundary  and  W [x,  —  o)  yt  0.  The  boundary  pressure  enters  in  a  term  similar  to  that 
for  any  other  particle  j. 

Does  this  expression  conserve  momentum?  If  we  multiply  by  m,  and  sum, 

we  find 


6 


d 

dt 


Y  m*u« 


PoJ2mt 


—  I  + 1  -j 

p 2 
+ 


/p\ 

Uv. 


p 

^2 


Y  Y1  m*mJ 


<  ) 


-  o) 

W(x,  -  X) 
dfE  (z,  —  Zj) 


dz{ 


(12) 


The  last  term  in  the  equation  vanishes  because  dW/dz,  is  antisymmetric  in  i  and 
j.  To  get  a  form  where  P,  and  P;  appear  symmetrically  was  the  main  point  in  the 
manipulation  of  P  and  p  in  Eq.  7.  This  insures  the  vanishing  of  the  last  term  in  Eq.  12 
which  is  necessary  for  conservation. 

Now  what  about  the  boundary  terms?  From  Eq.  3,  we  have 


-  o)  =  -p0  , 


(13) 


so  that  Eq.  12  becomes 


jt £ mtvi  =  \p0  +  poYmi(ji) 

-\Px-exY.mi\^W(xl-X).  (14) 

Here  P0  and  Px  are  the  given  boundary  values  of  the  external  pressure.  The  other 
terms  are  the  interpolated  values  of  the  pressure  at  the  boundary  given  by  the  internal 
calculation.  If  this  calculation  is  accurate  (apart  from  the  factor  of  1/2),  it  must  yield 


\po  -  PoYm'(j2  )  W ix'  ~  °)  ’ 


(15 


and 


7 


lPx^PxEm>(r)  W(xi-X), 

L  .  \P  ! i 


(16) 


giving  for  Eq.  12 


d_ 

dt 


E  rn,Vi  =  P0  -  Px  ■ 


(17) 


A  similar  analysis  is  possible  for  Monaghan’s  other  conservative  form  of  the 
momentum  equation.  First  write 


dv  _  l^P  _  2P1/Z  dP1'2 

dt  p  dx  p  dx 

Then  the  kernel  estimate  of  the  momentum  equation  becomes 


dvk 

dt 


Evaluating  the  integral  by  the  particle  method  gives 


(19) 


dvi 

dt 


p  1/2' 


-  2  E  m; 


[p}/2W(xt  -  o)  -  Plx/2W(Xi  -  X)] 
Pt1/2P;/2  dW(x,  ~  x,) 


PiPj 


dx, 


(20) 


where,  again,  we  note  the  fortuitous  appearance  of  the  factor  of  2  in  the  boundary 
terms  which  will  compensate  the  falloff  represented  by  Eq.  3. 
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After  multiplying  by  m  and  summing,  we  obtain 


Cl  r — \ 

—  >  m.u, 

dt  4- 


2 P'J2T.m<  [pf)  W'(x.-o) 

-2P'J2Y.m,  (— )  W(x<-X) 
i  \  P  J  i 


P!/2P}12  dW(x,  -x}) 


where  because  of  the  symmetry  in  t  and  j,  the  last  term  is  zero;  and  since 


W{Xi  -o), 

\pxl  -Hm.  w(Xi  ~  x )  - 


we  get  the  conservation  of  momentum, 


=  p.  -  PX  . 


It  is  curious  that  in  both  of  these  forms  the  boundary  terms  enter  in  such  a 
way  that  the  factor  of  1/2  from  the  normalization  anomaly  at  the  boundary  is  required 
for  a  conservative  momentum  equation.  This  is  not  true  of  all  forms,  however.  For 
example,  the  first  expression  for  the  pressure  gradient  in  Eq.  7  does  not  produce  a 
boundary  term  with  the  factor  of  2  present. 


SECTION  4 


EQUIVALENT  FINITE  DIFFERENCE  FORMS 


If  particles  are  assumed  equally  spaced,  it  is  possible  to  derive  the  equivalent 
finite  difference  form  which  the  kernel  estimate  represents.6  It  is  instructive  to  look  at 
these  forms,  particularly  in  the  neighborhood  of  the  boundary,  to  get  a  feeling  for  the 
accuracy  of  the  particle  method. 

For  the  purpose  of  illustration,  consider  the  kernel  shown  together  with  its 
derivative  in  Fig.  3, 2,6 


W  [x  —  x') 


1  (2 

*>  1  i  | 

7  -  - 

u2  +  -  u  3 

o  < 

h  \3 

2  1  / 

S<2- 

M)3 

i  < 

0 

2  < 

u  |  <  1 
u|  <  2 


where  u  =  (x  —  x')/h ,  and 


dW  _  1  dW 
dx  h  du 


(24) 


(25) 


Now,  consider  a  series  of  equally  spaced  equal  mass  particles  far  from  any 
boundary  as  shown  in  Fig.  A.  The  momentum  equation  in  the  form  of  Eq.  11  is 
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Figure  3.  The  smoothing  kernel  given  by  Eq.  24  and  its 
(continuous)  first  derivative. 
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since  dW,7  j dxx  is  zero  >.xcepi  for  j =  i  —  1  and  j  =  j  +  1  as  can  be  seen  from  Fig. 3. 
From  Eqs.  24-25  we  have 


dW (x,  -  I,)  _  I 
dx{ 


-1 

2 h? 

0 


j  =  i  ~  1 
j  =  * 


1 


j  =  i  +  l, 

which  leads  to  Eq.  26  above.  This  simplifies  to  the  central  difference  formula 


(2i 


dvi 

dt 


1 

2 h 


-) 


P 

,  p , 


»-  i  J 


where  we  have  used  m/p  —  h,  for  equally  spaced  particles. 


(28) 


Let  us  now  consider  the  case  where  the  boundary  falls  within  h  to  the  right  of 
particle  i.  In  this  case,  the  j  =  t  +  1  term  is  missing  but  the  boundary  term  is  present. 
Specifically,  let  us  assume  the  boundary  X  occurs  ~  0.6 h  to  the  right  of  particle  i  where 
W{ x,  -  X)  ~  1/2 h.  Eq.  11  then  becomes 


dv, 

dt 


12 


(29) 


!.•••••  i 

I  i-2  i-1  i  i+1  i+2  , 


(a) 


Figure  4(a).  A  series  of  equally-spaced  particles  far  from  a  boundary. 


(b) 


Figure  4(b).  Equivalent  difference  form  for  the  pressure  gradient  in  the 
neighborhood  of  the  boundary. 
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since  W  (x<  —  X)  =  1/2  h  and  the  j  =  »  + 1  term  is  missing.  Using  h  =  m/p,  this  reduces 
to 


dv{ 

dt 


(30) 


which  is  a  fairly  good  approximation.  As  shown  in  Fig.  4  it  gives  a  central  difference 
form  based  on  Px  situated  2h  from  P,_i,  although  it  was  assumed  that  the  boundary 
was  ~  1.6 h  from  x^_ j.  All  finite  difference  algorithms  are  uncertain  to  within  Ax  =  h, 
so  these  boundary  terms  are  probably  as  good  as  one  would  expect  of  most  difference 
methods. 
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SECTION  5 


WORK  TERMS  IN  THE  ENERGY  EQUATION 


We  now  consider  the  work  terms  which  appear  in  the  energy  equation  in  the 
presence  of  boundary  pressures.  In  order  to  insure  conservation,  this  term  must  be 
expressed  in  a  form  that  is  complimentary  to  the  right-hand  side  of  the  momentum 
equation  so  that  the  kinetic  energy  appears  correctly. 

We  write  the  work  terms  as 


dt 

dt 


Pdv 
p  dx 


d_ 

dx 


(31) 


where  c  is  the  internal  energy  per  unit  mass.  Taking  the  kernel  estimate  of  Eq.  31 
yields 


deK 

dt 


L 


+ 


*  A 

dx'  ^ 

x  JL 

dx' 


Pv 


W(x  —  x')dx' 


vk  ! 


—  \  W(x  —  x)dx  , 
v  P 


(32) 


where  the  second  term  has  been  linearized  as  before.  This  term  is  now  rewritten  using 
the  kernel  estimate  for  the  momentum  equation, 


dvK 

dt 


[x  d  (P\ 

Jo  dx1  \  p  ) 


W  (i  —  x)dx 


(33) 
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If  we  multiply  by  vK,  we  get 


which  gives  for  the  energy  equation 


(3-0 


deK 

dt 


(35) 


Next,  integrating  by  parts  we  obtain 


dtK  dvK 

- y  u,  — 

dt  K  dt 


(tH  -f( 


*  (Pv\  dw 


p  I  dx 


( x  -  x)dx' 


[pW]x_(?l\  txn™L 


j  J  P^{x-X')dX' 


(36) 


When  this  expression  is  evaluated  by  the  particle  method,  we  get  the  energy 
equation  expressed  in  terms  of  total  energy, 


df,  d  / 1  2\ 
dt  '  dt  V2  V 


PoW{Xi 


pxW(x,  - 


X) 


dW  (i,  -  Xj) 

dxt 


(37) 
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We  can  verify  that  Eq.  37  insures  conservation  if  we  multiply  by  m,  and  surn, 


d 

dt 


•f 


l._  ,  ( Pv\  W(xi-o) 

-(Pv)o  +  PoY.m>  I  —  ) - — - 

2  ,  V  P  )  i  Pi 

lln  S  V-  [Pv\W(x,-X ) 

-  ^{Pv)x  -  —  - - - 

2  V  P  /  Pi 


-  E  E  m<mj 


'  j 


P  V 


dWjXj  -  x j) 

dxt 


(38) 


The  left-hand  side  is  the  time  rate  of  change  of  total  energy,  internal  plus  kinetic.  The 
double  sum  on  the  right  is  zero  because  of  the  antisymmetry  of  dWtJ/dxt.  And  as 
before,  we  have 


-  2(Pv^°  ’ 


2(C»)x, 


(:») 


giving 


^  Em<(f,  f  ^,2)  =  (Pv)o  -  (Pv)x  ■  (10) 

»  ** 

'I'he  energy  Eq.  37  is  expressed  in  terms  of  the  total  energy.  After  the  particles 
are  moved  using  the  momentum  equation,  the  change  in  kinetic  energy  is  known.  Then 
Eq.  37  can  be  solved  f.r  the  internal  energy  In  some  cases,  it  may  be  desirable  to 
have  the  energy  equation  expressed  in  terms  of  the  internal  energy  alone. 
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The  explicit  appearance  of  the  kinetic  energy  term  can  he  eliminated  using 
the  momentum  Eq.  11, 


pcW  (x, 


°) 


P_ 


KP‘ 

- 


f  Pv' 


PxW(x,  -  X) 


Vi  + 


Pv\ 

p 2  J. 


dW{xi  -  x_,) 
dx, 


('ll) 


After  subtraction,  this  leads  to  a  form  of  the  energy  equation  where  the  kinetic  energy 
does  not  appear  explicitly, 


de, 

dt 


°) 


-  [ux  -  V, 


-  L 

J 


W{Xi-X) 

\pJx 


dW  (x,  —  x}) 


; 


dx. 


(•12) 


In  summary,  we  have  developed  a  conservative  expression,  Eq.  37,  for  the 
work  done  on  a  fluid  by  pressure  forces  that  is  complementary  to  the  momentum  Eq.  1 1 
in  the  presence  of  externally-applied  (boundary)  pressures.  These  two  equations  insure 
the  conservation  of  momentum  and  energy. 
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SECTION  6 


HEAT  CONDUCTION  TERMS  IN  THE  ENERGY 

EQUATION 


Let  us  now  consider  the  addition  of  thermal  heat  conduction  terms  in  the 
energy  equation.  This  will  be  especially  important  for  ICF  applications.  Neglecting 
other  terms  for  the  moment,  the  energy  equation  becomes 


d_i  =  ±±(k*£\ 

dt  p  d  x  \  dx  ) 

_1  dF 

p  dx 

where  K  is  the  thermal  conductivity,  T  the  temperature  and  F  the  heat  flux.  We 
shall  approach  this  equation  in  the  same  way  as  the  momentum  equation  and  obtain  a 
completely  parallel  formulation. 

Following  Section  3,  we  write 

dc  d  f  F\  F  dp 

dt  dx  \  p  )  p2  dx 

which  has  the  following  kernel  estimate 


(-13) 


d(  k 
dt 


After  integration  by  parts,  one  obtains 


(-15 


19 


deK 

dt 


F  1*  rx  l  F\  dW 

-W(x-x')  ~  -P-li-i'W 

.P  Jo  Jo  \p }  dxK  ’ 


F_ 


[pW(x  -  x')\0 


?U 


x  dW(x-x')^, 


dx 


dx'  , 


which  in  particle  form  becomes 


(•16) 


de{ 

dt 


+ 

0 


p0W  (x,  -  o) 


pxW{xx  -  X) 


Y. 


m. 


dW  (x,  —  Xj) 
dxx 


(47) 


It  is  clear  that  energy  conservation  follows  in  the  same  way  as  demonstrated  above. 

Thus,  in  Eq.  47  we  have  the  energy  equation  with  heat  conduction  in  a 
form  completely  analogous  to  the  momentum  equation,  Eq.  11.  The  boundary  values 
for  the  fluxes,  F0  and  Fx  appear  in  a  similar  way  to  the  boundary  pressures.  And 
since  experience  has  proven  Eq.  11  to  be  a  successful  particle  form  of  the  momentum 
equation,  one  would  expect  similar  results  for  Eq.  47. 

There  remains  only  to  calculate  the  fluxes  for  use  in  Eq.  47.  We  write  F  in 

the  form 


F  _  l  dT 
p  p  dx 


K 


K^dp 
p'dx  ’ 


dx  \p 

which,  by  a  now-familiar  procedure,  leads  to  the  particle  equation, 


(48) 
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F 

P 


(49) 


So  the  temperatures  define  the  fluxes  through  Eq.  49  and  the  fluxes  determine  the 
energy  transport  through  Eq.  47. 

It  is  possible  to  develop  an  equivalent  treatment  of  heat  conduction  based 
on  the  product  of  square  root  approach  used  in  Eq.  18,  but  this  is  unlikely  to  be  any 
better  than  the  method  already  given.  Other  approaches  have  been  proposed  for  heat 
conduction  (which  neglect  boundary  terms,  however).  The  method  of  Brookshaw  is 
discussed  in  the  following  section,  and  a  variation  of  Brookshaw’s  method  by  Monaghan 
gives  similar  results  in  test  calculations.5 
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SECTION  7 


CRITIQUE  OF  BROOKSHAW’S  METHOD  OF 
CALCULATING  THERMAL  DIFFUSION 


In  a  1985  paper,6  Brookshaw  discusses  several  approaches  to  adding  radiative 
heat  diffusion  to  SPII.  He  describes  Lucy’s  early  attempt7  at  solving  the  energy  equation 


pT 


dS_ 

dt 


F 


by  using  the  particle  method 


dF 

dx 

4  acT3  dT 
3Kp  dx 


(50) 


Fi 


( 4acT3\  y'(T]  dW(x<~xi) 

~  V  3 Kp  jr^rKp),  dxi 


(51) 


It  is  stated  that  since  W(x{  -  x,)  ~ >  0  at  the  boundary  of  the  fluid,  T  <■  0, 
and  one  has  an  insulating  boundary  condition.  The  statement  is  made,  .  .so  there  is 
no  reason  to  assume  that  the  SPH  equation  will  automatically  take  into  account  the 
correct  boundary  condition.” 


It  is  easy  to  see  where  Lucy  and  Brookshaw  go  wrong  here.  In  order  to  get  the 
particle  equations,  Eq.  51,  they  had  to  integrate  by  parts  and  drop  the  very  boundary 
terms  they  are  so  concerned  about.  For  example 


fX  Qf 

MW(I  ~  X')dz 

v  rx‘  dW 

-\FW(x-x')\0  -  Jo  F  —  [x-x')dx, 


(52) 
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which  becomes  in  particle  form 


F0W(Xi  -  o)  -  FxW(x,  -  X) 


(53) 


Even  though  boundary  terms  are  present  here,  this  is  not  a  good  form  to 
use,  since  it  does  not  guarantee  conservation.  It  is  better  to  use  the  method  developed 
in  Section  6.  Test  calculations  showed  Eq.  51  was  inaccurate  near  the  boundary  and 
also  the  method  was  sensitive  to  particle  positions.  Brookshaw  presented  a  different 
approach  to  calculating  heat  diffusion  based  on  a  difference  algorithm  that  he  claimed 
gives  much  smoother  results  than  Eq.  51.  The  method  given  in  Section  6  should  be 
accurate  near  the  boundaries,  but  it  may  be  that  Brookshaw’s  approach  would  give 
added  smoothness.  For  this  reason,  we  shall  attempt  to  derive  Brookshaw’s  algorithm 
with  boundary  terms  included. 

Write  the  energy  equation  as 


de 


P 


dt 


and  form  the  kernel  estimate 


Pk 


deK 

dt 


rx  d_ 

o  dx' 


W{  x 


x')dx'  , 


(54) 


(55) 


where  curly  brackets  are  used  to  indicate  a  special  treatment  for  that  particular  term. 
Integrating  by  parts  yields 


Pk 


dju 

dt 


W  (x  -  x') 


+ 


[*(  dT\dW 

L  \Ksi’\^[x~z)dx 


(56) 
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Now,  Brookshaw  appeals  to  a  Taylor  series  expansion  to  write  a  difference 
form  for  the  term  in  curly  brackets. 

Since 


K(x')  =  K(x)-(x-x,)^  +  ... 

T(x')  =  r(x)-(x-x')g  +  i(x-x')! 


d2T 
dx 2 


+  ...  , 


(57) 


we  can  write 


tffj  }’  =  [*(*)  +  fi-(x')!  TMx  - 


=  2K{x)~-  -  (x  -  x') 

OX 


d2T  dKdT 

^  ax2  +  ax  ax 


+ 


(58) 


The  kernel  estimate,  Eq.  56,  becomes 


+  f  |JT(x)  +  irwi  r(l|  I  -  xVx’  .  (50) 

where  Brookshaw  only  discusses  the  integral  and  assumes  everything  goes  to  zero  at 
the  boundaries. 

It  is  interesting  to  see  something  like  2  K  under  the  integral  sign.  To  see 
how  the  difference  form  gives  the  intended  heat  conduction  term,  substitute  Eq.  58 
into  Eq.  59 


Pk 


du 

dt 


K 


a r 

dx 


hx  \  dx 


+  ...\W(x-x') 

a 


+ . . . 


dx 


\V{x 


x)dx  .  (60) 
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When  the  integral  is  performed,  term-by-term  cancellation  gives 


Pk 


deK 

dt 


+  0(h 2)  , 


(61) 


which  is  the  result  obtained  by  Brookshaw.  Including  the  boundary  terms  in  the 
integration  by  parts  does  not  change  the  result  given  by  series  expansion. 

Brookshaw’s  scheme  with  boundary  terms  can  now  be  written  as 


=  [<*(*)  +  g(x')}  r(*]  _  ^XV(s  -  *')] 


nix 
0 


+  L  +  K(x')}  T^x~_  x  ^  ~  x')dx'  ’ 


(62) 


which  gives  in  particle  form 


dt, 

dt 


--[Ki  +  K0 

Pi 


Ti  -  Tc 


T 

i 


X: 


W(Xi-o)  +  ~[K,+Kx  1 
o  Pi  x,  —  A 


W(Xi  -  X) 


+  mL - (x«  -  xj)  ■ 

j  PiPj  •£*  uXx 


(63) 


It  is  clear  that  Eq.  63  gives  conservation  since  the  summation  is  antisymmetric 
in  i  and  j  and  the  necessary  factors  of  2  are  present  in  the  boundary  terms.  Whether 
this  approach  to  heat  conduction  is  better  than  that  represented  by  Eqs.  47-49  can  only 
be  determined  by  numerical  testing.  Eq.  63  certainly  solves  the  problem  of  an  insulating 
boundary  condition,  since  any  particle  within  2 h  of  the  boundary  (which  presumably 
is  on  the  order  of  one  diffusion  length)  will  exchange  energy  with  the  environment 
icprcscnt.ed  by  the  temperatures  T0  and  Tx. 
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SECTION  8 


CRITIQUE  OF  SOME  MANIPULATIONS  IN  THE 

LITERATURE 


In  the  foregoing,  it  has  been  shown  that  boundary  terms  neglected  in  other 
derivations  of  SPH  algorithms2-7  are  essential  for  many  problems  and  particularly 
for  energy  transport  problems.  The  question  arises  as  to  how  these  terms  were  so 
consistently  neglected  over  the  last  decade. 

In  G ingold  and  Monaghan’s  paper,4  they  approach  the  particle  method  as 

follows: 


Kernel  Estimate  PK(x) 
Particle  Approx.  Pt 
dP 

The  Gradient  — — 
dxi 


J  P{x')W{x  —  x)dx  , 

J2  mi  (  ~  )  w (z*  ~  xj)  - 
j  J  j 


(64) 


However,  if  one  calculates  the  pressure  gradient  directly  from  the  kernel  estimate  one 
finds: 


Kernel  Estimate 
Integration 
Particle  Approx. 


dP 

dxi 


r*  dP 

=  l  a?w(x'x')dx'' 

v  fX  dW 

=  l^|f  +  /o  P{x')  —  (x-x’)dx', 


=  PxW{xi-X)-P0W{xi-o) 


+ 


j 


dW  (xi  —  x}) 
dxi 


(65) 


26 


This  shows  clearly  that  one  must  be  careful  of  taking  derivatives  through 
particle  approximations  to  kernel  estimates.  In  Monaghan’s  1982  paper,3  he  mentions 
integration  by  parts  but  says  he  is  assuming  W ,  the  function,  or  both  go  to  zero  on 
the  boundary.  This  is  better  than  what  one  sees  in  References  4  and  6,  because  it  is 
mathematically  defined  even  though  unnecessarily  restrictive. 

The  point  to  be  made  here  is  that  in  dealing  with  SPH  equations,  it  is  safer 
to  carry  everything  through  to  final  form  as  an  integral  equation  and  only  then  replace 
the  integrals  with  sums  over  particles. 
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SECTION  9 


EXTENSION  TO  HIGHER  DIMENSIONS 


Consider  first  the  kernel  estimate  for  the  momentum  equation,  Eq.  7.  In 
vector  notation,  this  becomes 


-  ^  Jv  V'pW(\f-  r'|,A)dV,  (60) 

where  we  assume  the  kernel  is  spherically  summetric,  extending  over  a  range  R  ~  2h 
about  the  point  r.  The  volume  of  integration  V  is  therefore  a  sphere  of  radius  U 
centered  on  fas  in  Fig.  5.  An  exception  to  this  is  when  f  is  within  range  of  a  boundary 
as  in  Fig.  6.  Then  integration  by  parts  yields  a  surface  term,  where  the  boundary 
pressure  at  the  surface  influences  the  particle  at  point  r. 

Integrating  by  parts  in  Eq.  66  gives 

-(£)  JBpW{\f~f\h)dS,~  ^  JvpVW(\r-r'\,h)dV ,  (67) 


where  the  surface  integrals  are  over  that  area  of  the  boundary  intercepted  by  t  he  non¬ 
zero  range  of  W .  The  particle  approximation  to  the  kernel  estimate  becomes 
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o 


Figure  5.  Range  of  influence  of  kernel  VF(|r  —  r'|)  about  the 
point  r. 


Figure  6.  Range  of  influence  of  W  intercepts  boundary 
plane  B  over  area  AS. 
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Figure  7.  The  area  on  a  plane  boundary  intercepted  by  the 
sphere  is  a  disk  of  radius  a. 

-Z>>  (^)  +(^)  VMK(|r;  -  r,| ,h),  (68) 

where  pB  and  PB  are  the  density  and  pressure  at  the  boundary  averaged  over  AS,,  rB 
is  the  normal  point  at  the  boundary  closest  to  r,  and  h  is  the  unit  vector  normal  to 
the  boundary  at  AS,,  all  as  shown  in  Fig.  6. 

The  calculation  of  AS,  is  fairly  easy  if  one  neglects  curvature  of  the  boundary, 
which  is  justified  in  most  cases  since  h  is  supposed  to  be  small.  In  three  dimensions, 
the  area  intercepted  by  a  sphere  of  radius  R  located  a  distance  d  =  |r  -  rfl|  along  the 
normal  to  a  plane  boundary,  as  illustrated  in  Fig.  7,  is 

ASj  =  7r  (f?2  —  |r  —  rB|2)  .  (69) 

In  2D  rectangular  geometry  (x,  y)  one  considers  a  unit  length  in  2.  Hence,  the  range 
of  influence  of  W  is  a  cylinder  of  radius  R  and  unit  length.  In  this  case,  which  is 
illustrated  in  Fig.  8,  the  area  intercepted  by  W  on  a  plane  boundary  is 

AS,  =  2(f?2-|r-rB!2),/2.  (70) 


30 


Figure  8.  The  area  on  a  plane  boundary  intercepted  by  the 
cylinder  of  unit  length  is  a  rectangle  of  height  2 a. 

This  formula  does  not  apply  to  the  situation  in  2D  cylindrical  geometry  where,  be¬ 
cause  of  the  symmetry,  a  particle  is  toroidal  in  shape.  There  does  not  appear  to  be 
a  simple  expression  that  covers  all  boundary  configurations  in  this  geometry.  In  the 
one-dimensional  cases  discussed  earlier,  one  considers  a  unit  length  in  each  of  the  two 
transverse  directions  giving  AS,  =  1.  Eq.  68  then  reduces  to  the  ID  expression  derived 
earlier.  The  area  AS,  in  these  three  cases  can  be  written  generally  as 

AS,  -  C{a)(R2  -  I*  -  rB|2)(“~1)/2  ,  (71) 

where  a  --  1,2,3  is  the  geometrical  order,  and  C(a)  =  (1,2, 7r),  respectively.  It  should 
be  noted  that  an  exception  to  Eq.  71  may  occur  at  the  corner  of  intersecting  boundaries. 

The  normalization  anomaly  discussed  in  Section  2  is  not  altered  in  higher 
geometries.  In  this  case,  one  has 


lim  [  W{  1 

r->rD  Jv 


r-r'i  ,fc)dV  =  -, 


(71 


and  for  a  uniform  distribution  of  particles  in  the  neighborhood  of  the  boundary, 


J2m,W(\fD  -  f,\,h)  =  ^ pc  . 

J 


(73) 
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Likewise,  the  boundary  value  of  the  pressure  or  any  similar  function  in  the  particle 
approximation  would  be  expressed  as 

-«.*)•  (74) 

In  order  to  demonstrate  conservation  of  momentum  from  Eq.  68,  one  obtains 
the  total  component  of  momentum  in  any  direction  x, 


»  — ♦ 

=  -pB52rni(PBASi-x)W(\ri-fB\,h) 

i  t 

^A$-xj  VV(|»\ -fB|,/i)  ,  (75) 

where  the  double  sum  vanishes  because  of  the  antisymmetry  in  i  and  j.  The  first  term 
is  the  particle  approximation  to  one-half  of  the  boundary  value  of  the  net  force  on  the 
system  of  particles, 


i  B PAS,  ■  i)„  ,  (76) 

and  by  the  same  arguments  presented  earlier  in  Section  3,  the  second  term  should  be 
a  good  approximation  to  this  expression  also.  Conservation  is  then  expressed  as 


E 


z=Y,{P&Siz)B. 

i 


(77) 
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By  a  similar  series  of  steps  the  work  terms  in  the  energy  equation  of  Section  5 
are  written 


dt,  d 
dt  n  dt 


AS,pBW(\r,  -  rB\,h) 


E 


m, 


VW  (jrq  —  r-j!,  h), 


(78) 


and  conservation  of  energy  follows  in  a  like  manner, 


Jt  Em*(«.  +  \vi)  =  ■  *>) b- 


(79) 


The  heat  flux  terms  in  Section  6  are  easily  generalized  to  the  form. 


dt, 

~dt 


•  A S,pBW (|f,  -  rB\,h) 


+ 

/ 


•  v^((f;  -fB\,h) , 


(80) 


with  the  fluxes  given  by 


A S,pBW{\r,  -  rB\,h) 


VW(\f,  -fjUh). 


(81) 
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SECTION  10 


SUMMARY  OF  NEW  SPH  ALGORITHMS 


A  suggested  approach  to  SPH  is  summarized  here  that  allows  boundary  values 
for  pressure,  temperature  and  heat  flux  to  be  specified.  This  allows  problems  to  be 
solved  where  the  fluid  is  driven  by  an  externally-applied  pressure,  temperature  or  both. 
It  also  allows  the  fluid  to  lose  energy  to  a  cooler  environment.  Although  the  diffusion 
treatment  is  derived  for  thermal  conduction,  radiation  which  is  in  equilibrium  with  the 
material  temperature  can  be  handled  in  the  same  way. 

The  momentum  equation  for  particle  t  is  written 


dvt 

dt 


pBW(\ft-fB\,k)ASt 


//,  +  \P\ 


VW(| Ti-f^h)  . 


(82) 


In  this  equation  PB  is  the  external  scalar  pressure  imposed  on  the  fluid,  specified  as 
problem  input.  The  quantity  pB  is  the  boundary  value  of  the  density  given  by 


PB  =  ^^2mjW{\ri-fB\,h).  (83) 

i 

The  factor  of  2  appears  as  in  Eq.  74  because  of  the  normalization  anomaly  at  the 
boundary. 


From  Eq.  82  we  see  that  any  particle  i  within  range  of  the  boundary,  i.e., 
|ri  ~  F/j|<  2h,  begins  to  be  influenced  by  the  external  pressure.  This  pressure  term 
enters  in  a  form  similar  to  that  of  any  neighboring  particle  except  W  appears  in  the 
boundary  term  rather  than  VVF.  This  equation  is,  therefore,  different  than  if  one 
assumed  that  a  fixed  particle  represented  the  boundary  with  a  pressure  defined  there. 
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The  work  term  which  appears  in  the  energy  equation  is  written 


d 

dt  dt 


(84) 


After  the  particles  are  moved  using  Eq.  82,  the  change  in  kinetic  energy  is  known  and 
Eq.  84  is  used  to  obtain  the  change  in  internal  energy.  From  the  new  values  of  p,  and 
t,,  the  temperatures  and  pressures  can  be  calculated. 

If  thermal  or  equilibrium  radiation  diffusion  is  present,  the  following  heat  flux 
terms  must  be  added  to  the  work  terms  in  the  energy  equation 


~  --  Work  Terms 

dt 


(85) 


The  heat  fluxes  Ft  are  defined  at  each  particle  position  just  as  the  pressures  and  tem¬ 
peratures  are.  The  fluxes  are  obtained  from  the  temperatures  by 


i 


(86) 


It  is  not  clear  how  best  to  implement  the  diffusion  equation  in  the  finite  dif¬ 
ference  approximation.  Some  numerical  experimentation  may  be  necessary.  However, 
the  above  set  of  SPH  equations  allows  for  boundary  conditions  on  P,  F  and  T  and 
conserves  energy,  momentum,  and  mass. 
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